Project: 1st Order Spatial Point Patterns Analysis Methods

pacman::p_load(readxl, httr, jsonlite, maptools, sf, raster, spatstat, spNetwork, rgdal, sp, tmap, tidyverse)

Geospatial Data Wrangling

For 1st order SPPA: airbnb_sf, mpsz_sf, sg_sf For 2nd order SPPA: same as above For NetSPPA: roadnetwork_sf For K function, K cross function, LCLQ: sevenele_sf, busstop_sf, hotel_sf, mrt_sf, mall_sf, attr_sf, unis_sf

read the 6 shp, kml, geojson files

busstop_sf <- st_read(dsn = "data/geospatial/busstop-022023-shp", layer="BusStop")

Reading layer `BusStop' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\busstop-022023-shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 5159 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21

mpsz_sf <- st_read(dsn = "data/geospatial/MPSZ-2019-shp", layer="MPSZ-2019")

Reading layer `MPSZ-2019' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\MPSZ-2019-shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS:  WGS 84

network_sf <- st_read(dsn = "data/geospatial/roadnetwork-022023-shp", layer="RoadSectionLine")

Reading layer `RoadSectionLine' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\roadnetwork-022023-shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 15300 features and 2 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 3248.394 ymin: 22755.1 xmax: 50043.6 ymax: 50135.95
Projected CRS: SVY21

attr_sf <- st_read(dsn = "data/geospatial/touristattractions-2017-shp", layer="TOURISM")

Reading layer `TOURISM' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\touristattractions-2017-shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 107 features and 19 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -1.797693e+308 ymin: -1.797693e+308 xmax: 43659.54 ymax: 47596.73
Projected CRS: SVY21

mrt_sf <- st_read(dsn = "data/geospatial/mrt-112022-shp", layer="Train_Station_Exit_Layer") 

Reading layer `Train_Station_Exit_Layer' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\mrt-112022-shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 562 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 6134.086 ymin: 27499.7 xmax: 45356.36 ymax: 47865.92
Projected CRS: SVY21

hotel_sf <- st_read("data/geospatial/hotels-2021-kml/hotel-locations.kml")

Reading layer `HOTELS' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\hotels-2021-kml\hotel-locations.kml' 
  using driver `KML'
Simple feature collection with 422 features and 2 fields
Geometry type: POINT
Dimension:     XYZ
Bounding box:  xmin: 103.6351 ymin: 1.245797 xmax: 103.9891 ymax: 1.419278
z_range:       zmin: 0 zmax: 0
Geodetic CRS:  WGS 84

sg_sf <- st_read(dsn = "data/geospatial/costal-shp", layer="CostalOutline")

Reading layer `CostalOutline' from data source 
  `C:\valtyl\IS415-GAA\Project\data\geospatial\costal-shp' using driver `ESRI Shapefile'
Simple feature collection with 60 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2663.926 ymin: 16357.98 xmax: 56047.79 ymax: 50244.03
Projected CRS: SVY21

read mall csv and convert to sf

mall_csv <- read_csv("data/geospatial/shoppingmalls-2019-csv/mall_coordinates_updated.csv")
mall_sf <- st_as_sf(mall_csv, coords = c("longitude", "latitude"), crs=4326)

mall_sf

Simple feature collection with 184 features and 2 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 103.6784 ymin: 1.263797 xmax: 103.9897 ymax: 1.448227
Geodetic CRS:  WGS 84
# A tibble: 184 × 3
    ...1 name                               geometry
 * <dbl> <chr>                           <POINT [°]>
 1     0 100 AM                  (103.8435 1.274588)
 2     1 112 KATONG              (103.9051 1.305087)
 3     2 313@SOMERSET            (103.8377 1.301385)
 4     3 321 CLEMENTI             (103.765 1.312025)
 5     4 600 @ TOA PAYOH          (103.851 1.334042)
 6     5 888 PLAZA               (103.7953 1.437131)
 7     6 AMK HUB                 (103.8485 1.369223)
 8     7 ADMIRALTY PLACE         (103.8018 1.439881)
 9     8 ALEXANDRA RETAIL CENTRE (103.8014 1.273843)
10     9 ANCHORPOINT             (103.8056 1.288935)
# … with 174 more rows

read airbnb csv and convert to sf

airbnb_csv <- read_csv("data/geospatial/airbnb-2022-csv/listings.csv")
airbnb_sf <- st_as_sf(airbnb_csv, coords = c("longitude", "latitude"), crs=4326)

get_coords <- function(add_list){
  
  # Create a data frame to store all retrieved coordinates
  postal_coords <- data.frame()
  
  # loop to go through each address in the list  
  for (i in add_list){
    #print(i)
    
    # response from OneMap API
    r <- GET('https://developers.onemap.sg/commonapi/search?',
           query=list(searchVal=i,
                     returnGeom='Y',
                     getAddrDetails='Y'))
    data <- fromJSON(rawToChar(r$content))
    found <- data$found
    res <- data$results
    
    # Create a new data frame for each address
    new_row <- data.frame()
    
    # If single result, append 
    if (found == 1){
      postal <- res$POSTAL 
      lat <- res$LATITUDE
      lng <- res$LONGITUDE
      new_row <- data.frame(address= i, postal = postal, latitude = lat, longitude = lng)
    }
    
    # If multiple results, drop NIL and append top 1
    else if (found > 1){
      # Remove those with NIL as postal
      res_sub <- res[res$POSTAL != "NIL", ]
      
      # Set as NA first if no Postal
      if (nrow(res_sub) == 0) {
          new_row <- data.frame(address= i, postal = NA, latitude = NA, longitude = NA)
      }
      
      else{
        top1 <- head(res_sub, n = 1)
        postal <- top1$POSTAL 
        lat <- top1$LATITUDE
        lng <- top1$LONGITUDE
        new_row <- data.frame(address= i, postal = postal, latitude = lat, longitude = lng)
      }
    }

    else {
      new_row <- data.frame(address= i, postal = NA, latitude = NA, longitude = NA)
    }
    
    # Add the row to the main dataframe
    postal_coords <- rbind(postal_coords, new_row)
  }
  return(postal_coords)
}

read 7-11 xlsx and get coords

Show the code!

sevenele_xlsx <- read_xlsx("data/geospatial/711stores-032023-xlsx/7-11_convenience_stores.xlsx")

sevenele_list <- sort(unique(sevenele_xlsx$Postal))
sevenele_coords <- get_coords(sevenele_list)
sevenele_coords[(is.na(sevenele_coords$postal) | is.na(sevenele_coords$latitude) | is.na(sevenele_coords$longitude) | sevenele_coords$postal=="NIL"), ]

fix 711 coords

Show the code!

sevenele_coords[64,]$postal <- "098975"
sevenele_coords[64,]$latitude <- "1.255169"
sevenele_coords[64,]$longitude <- "103.818508"

sevenele_coords[455,]$postal <- "828694"
sevenele_coords[455,]$latitude <- "1.420668"
sevenele_coords[455,]$longitude <- "103.912257"

check if fixed

Show the code!

sevenele_coords[(is.na(sevenele_coords$postal) | is.na(sevenele_coords$latitude) | is.na(sevenele_coords$longitude) | sevenele_coords$postal=="NIL"), ]

get 711 rds

Show the code!

sevenele_df <- left_join(sevenele_xlsx, sevenele_coords, by=c('Postal' = 'address'))

sevenele_rds <- write_rds(sevenele_df, "data/geospatial/711stores-032023-xlsx/sevenele_rds.rds")

get 711 sf

seveneles <- read_rds("data/geospatial/711stores-032023-xlsx/sevenele_rds.rds")
sevenele_sf <- st_as_sf(seveneles,
                    coords = c("longitude", 
                               "latitude"),
                    crs=4326) %>%
  st_transform(crs = 3414)

read universities and colleges xlsx and get coords

Show the code!

unis_xlsx <- read_xlsx("data/geospatial/universitiescolleges-2023-xlsx/universities_and_colleges.xlsx")

unis_list <- sort(unique(unis_xlsx$name))
unis_coords <- get_coords(unis_list)
unis_coords[(is.na(unis_coords$postal) | is.na(unis_coords$latitude) | is.na(unis_coords$longitude) | unis_coords$postal=="NIL"), ]

fix uni coords

Show the code!

unis_coords[1,]$postal <- "179104"
unis_coords[1,]$latitude <- "1.29006"
unis_coords[1,]$longitude <- "103.8497"

unis_coords[2,]$postal <- "229469"
unis_coords[2,]$latitude <- "1.31648392303"
unis_coords[2,]$longitude <- "103.873964955"

unis_coords[4,]$postal <- "408601"
unis_coords[4,]$latitude <- "1.3196562"
unis_coords[4,]$longitude <- "103.8923565"

unis_coords[5,]$postal <- "059405"
unis_coords[5,]$latitude <- "1.28799"
unis_coords[5,]$longitude <- "103.84692"

unis_coords[6,]$postal <- "188946"
unis_coords[6,]$latitude <- "1.3462717"
unis_coords[6,]$longitude <- "103.95302009"

unis_coords[7,]$postal <- "329162"
unis_coords[7,]$latitude <- "1.3270276324741"
unis_coords[7,]$longitude <- "103.85336279869"

unis_coords[10,]$postal <- "188647"
unis_coords[10,]$latitude <- "1.2990616"
unis_coords[10,]$longitude <- "103.8494881"

unis_coords[11,]$postal <- "088383"
unis_coords[11,]$latitude <- "1.28013"
unis_coords[11,]$longitude <- "103.84193"

unis_coords[12,]$postal <- "049319"
unis_coords[12,]$latitude <- "1.28469"
unis_coords[12,]$longitude <- "103.85269"

unis_coords[14,]$postal <- "138676"
unis_coords[14,]$latitude <- "1.299122"
unis_coords[14,]$longitude <- "103.788025"

unis_coords[16,]$postal <- "228095"
unis_coords[16,]$latitude <- "1.30209070954"
unis_coords[16,]$longitude <- "103.849560613"

unis_coords[17,]$postal <- "188306"
unis_coords[17,]$latitude <- "1.3001799"
unis_coords[17,]$longitude <- "103.84924"

unis_coords[19,]$postal <- "069542"
unis_coords[19,]$latitude <- "1.2750794"
unis_coords[19,]$longitude <- "103.8462627"

unis_coords[20,]$postal <- "148951"
unis_coords[20,]$latitude <- "1.2970726"
unis_coords[20,]$longitude <- "103.8011673"

unis_coords[22,]$postal <- "189655"
unis_coords[22,]$latitude <- "1.30020700501"
unis_coords[22,]$longitude <- "103.851710836"

unis_coords[26,]$postal <- "248922"
unis_coords[26,]$latitude <- "1.31532"
unis_coords[26,]$longitude <- "103.77994"

unis_coords[30,]$postal <- "599491"
unis_coords[30,]$latitude <- "1.326"
unis_coords[30,]$longitude <- "103.79999"

unis_coords[37,]$postal <- "550266"
unis_coords[37,]$latitude <- "1.35324"
unis_coords[37,]$longitude <- "103.87145"

unis_coords[38,]$postal <- "139660"
unis_coords[38,]$latitude <- "1.3079"
unis_coords[38,]$longitude <- "103.777"

check fix

Show the code!

unis_coords[(is.na(unis_coords$postal) | is.na(unis_coords$latitude) | is.na(unis_coords$longitude) | unis_coords$postal=="NIL"), ]

get uni rds

Show the code!

unis_df <- left_join(unis_xlsx, unis_coords, by=c('name' = 'address'))
unis_rds <- write_rds(unis_df, "data/geospatial/universitiescolleges-2023-xlsx/unis_rds.rds")

get unis sf

unis <- read_rds("data/geospatial/universitiescolleges-2023-xlsx/unis_rds.rds")
unis_sf <- st_as_sf(unis,
                    coords = c("longitude", 
                               "latitude"),
                    crs=4326) %>%
  st_transform(crs = 3414)

remove unnecessary columns (except mpsz_sf, network_sf)

attr_sf <- attr_sf %>% select(c(6))
busstop_sf <- busstop_sf %>% select(c(1))
hotel_sf <- hotel_sf %>% select(c(1))
mall_sf <- mall_sf %>% select(c(2))
mrt_sf <- mrt_sf %>% select(c(1))
sevenele_sf <- sevenele_sf %>% select(c(1))
unis_sf <- unis_sf %>% select(c(1))
airbnb_sf <- airbnb_sf %>% select(c(2,6,7,8))

check for invalid geometries

length(which(st_is_valid(attr_sf) == FALSE))

[1] 0

length(which(st_is_valid(busstop_sf) == FALSE))

[1] 0

length(which(st_is_valid(hotel_sf) == FALSE))

[1] 0

length(which(st_is_valid(mall_sf) == FALSE))

[1] 0

length(which(st_is_valid(mrt_sf) == FALSE))

[1] 0

length(which(st_is_valid(sevenele_sf) == FALSE))

[1] 0

length(which(st_is_valid(unis_sf) == FALSE))

[1] 0

length(which(st_is_valid(mpsz_sf) == FALSE))

[1] 6

length(which(st_is_valid(network_sf) == FALSE))

[1] 0

length(which(st_is_valid(airbnb_sf) == FALSE))

[1] 0

length(which(st_is_valid(sg_sf) == FALSE))

[1] 1

fix mpsz invalid geometry and check again

mpsz_sf <- st_make_valid(mpsz_sf)
length(which(st_is_valid(mpsz_sf) == FALSE))

[1] 0

fix sg invalid geometry and check again

sg_sf <- st_make_valid(sg_sf)
length(which(st_is_valid(sg_sf) == FALSE))

[1] 0

dealing w missing values

attr_sf[rowSums(is.na(attr_sf))!=0,]

Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] PAGETITLE geometry 
<0 rows> (or 0-length row.names)

busstop_sf[rowSums(is.na(busstop_sf))!=0,]

Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] BUS_STOP_N geometry  
<0 rows> (or 0-length row.names)

hotel_sf[rowSums(is.na(hotel_sf))!=0,]

Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
[1] Name     geometry
<0 rows> (or 0-length row.names)

mall_sf[rowSums(is.na(mall_sf))!=0,]

Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
# A tibble: 0 × 2
# … with 2 variables: name <chr>, geometry <GEOMETRY [°]>

mrt_sf[rowSums(is.na(mrt_sf))!=0,]

Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] stn_name geometry
<0 rows> (or 0-length row.names)

sevenele_sf[rowSums(is.na(sevenele_sf))!=0,]

Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21 / Singapore TM
# A tibble: 0 × 2
# … with 2 variables: Name <chr>, geometry <GEOMETRY [m]>

unis_sf[rowSums(is.na(unis_sf))!=0,]

Simple feature collection with 0 features and 1 field
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21 / Singapore TM
# A tibble: 0 × 2
# … with 2 variables: name <chr>, geometry <GEOMETRY [m]>

mpsz_sf[rowSums(is.na(mpsz_sf))!=0,]

Simple feature collection with 0 features and 6 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
[1] SUBZONE_N  SUBZONE_C  PLN_AREA_N PLN_AREA_C REGION_N   REGION_C   geometry  
<0 rows> (or 0-length row.names)

network_sf[rowSums(is.na(network_sf))!=0,]

Simple feature collection with 15300 features and 2 fields
Geometry type: MULTILINESTRING
Dimension:     XY
Bounding box:  xmin: 3248.394 ymin: 22755.1 xmax: 50043.6 ymax: 50135.95
Projected CRS: SVY21
First 10 features:
   RD_CD            RD_CD_DESC                       geometry
1   <NA>       SEA BREEZE ROAD MULTILINESTRING ((41260.55 ...
2   <NA>       HILLVIEW AVENUE MULTILINESTRING ((20404.4 3...
3   <NA>     SERANGOON CENTRAL MULTILINESTRING ((32508.6 3...
4   <NA> CHANGI SOUTH STREET 3 MULTILINESTRING ((42949.75 ...
5   <NA>           PANDAN ROAD MULTILINESTRING ((18683.23 ...
6   <NA>      WEST COAST GROVE MULTILINESTRING ((19396.72 ...
7   <NA>          NEYTHAL ROAD MULTILINESTRING ((12917.95 ...
8   <NA>         WILLOW AVENUE MULTILINESTRING ((32638.54 ...
9   <NA> CHOA CHU KANG NORTH 5 MULTILINESTRING ((18182.78 ...
10  <NA>       CASSIA CRESCENT MULTILINESTRING ((33711.34 ...

airbnb_sf[rowSums(is.na(airbnb_sf))!=0,]

Simple feature collection with 0 features and 4 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS:  WGS 84
# A tibble: 0 × 5
# … with 5 variables: name <chr>, neighbourhood <chr>, room_type <chr>,
#   price <dbl>, geometry <GEOMETRY [°]>

sg_sf[rowSums(is.na(sg_sf))!=0,]

Simple feature collection with 0 features and 4 fields
Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
Projected CRS: SVY21
[1] GDO_GID    MSLINK     MAPID      COSTAL_NAM geometry  
<0 rows> (or 0-length row.names)

verify coordinate system

st_crs(attr_sf)

Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

st_crs(busstop_sf)

Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

st_crs(hotel_sf)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

st_crs(mall_sf)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

st_crs(mrt_sf)

Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

st_crs(sevenele_sf)

Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

st_crs(unis_sf)

Coordinate Reference System:
  User input: EPSG:3414 
  wkt:
PROJCRS["SVY21 / Singapore TM",
    BASEGEOGCRS["SVY21",
        DATUM["SVY21",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4757]],
    CONVERSION["Singapore Transverse Mercator",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["northing (N)",north,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["easting (E)",east,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["Cadastre, engineering survey, topographic mapping."],
        AREA["Singapore - onshore and offshore."],
        BBOX[1.13,103.59,1.47,104.07]],
    ID["EPSG",3414]]

st_crs(mpsz_sf)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

st_crs(network_sf)

Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

st_crs(airbnb_sf)

Coordinate Reference System:
  User input: EPSG:4326 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

st_crs(sg_sf)

Coordinate Reference System:
  User input: SVY21 
  wkt:
PROJCRS["SVY21",
    BASEGEOGCRS["SVY21[WGS84]",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Transverse Mercator",
            ID["EPSG",9807]],
        PARAMETER["Latitude of natural origin",1.36666666666667,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",103.833333333333,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",1,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",28001.642,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",38744.572,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

transform coordinate system

# with st_set_crs(), we can assign the appropriate ESPG Code
attr_sf <- st_set_crs(attr_sf, 3414)
busstop_sf <- st_set_crs(busstop_sf, 3414)
mrt_sf <- st_set_crs(mrt_sf, 3414)
network_sf <- st_set_crs(network_sf, 3414)
sg_sf <- st_set_crs(sg_sf, 3414)

# with st_transform(), we can change from one CRS to another
mpsz_sf <- st_transform(mpsz_sf, crs=3414)
hotel_sf <- st_transform(hotel_sf, crs=3414)
mall_sf <- st_transform(mall_sf, crs=3414)
airbnb_sf <- st_transform(airbnb_sf, crs=3414)

# sevenele_sf and unis_sf are in the correct CRS and ESPG code

Preparing for 1st Order SPPA

convert sf to sp’s spatial class

airbnb <- as_Spatial(airbnb_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(sg_sf)

convert sp’s spatial into generic sp format

airbnb_sp <- as(airbnb, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")

airbnb

class       : SpatialPointsDataFrame 
features    : 3037 
extent      : 7406.989, 44064, 25352.92, 48321.84  (xmin, xmax, ymin, ymax)
crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs 
variables   : 4
names       :                                               name, neighbourhood,       room_type,  price 
min values  :      -Master bedroom with attached bathroom at CBD,    Ang Mo Kio, Entire home/apt,      0 
max values  : 阿裕尼公寓普通房急需女搭房一名联系人:郑小姐 电话：,        Yishun,     Shared room, 135140 

convert generic sp format into spatstat’s ppp format

airbnb_ppp <- as(airbnb_sp, "ppp")
airbnb_ppp

Planar point pattern: 3037 points
window: rectangle = [7406.99, 44064] x [25352.92, 48321.84] units

check number of duplicate point events

sum(multiplicity(airbnb_ppp) > 1)

[1] 211

handle duplicate points:

airbnb_ppp_jit <- rjitter(airbnb_ppp, 
                             retry=TRUE, 
                             nsim=1, 
                             drop=TRUE)

check duplicated points:

any(duplicated(airbnb_ppp_jit))

[1] FALSE

create owin

sg_owin <- as(sg_sp, "owin")

combine airbnb point events obj and owin obj

airbnbSG_ppp = airbnb_ppp[sg_owin]

summary

summary(airbnbSG_ppp)

Planar point pattern:  3035 points
Average intensity 4.05347e-06 points per square unit

*Pattern contains duplicated points*

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Window: polygonal boundary
61 separate polygons (1 hole)
                  vertices         area relative.area
polygon 1               38  1.56140e+04      2.09e-05
polygon 2              735  4.69093e+06      6.27e-03
polygon 3               49  1.66986e+04      2.23e-05
polygon 4               76  3.12332e+05      4.17e-04
polygon 5             5141  6.36179e+08      8.50e-01
polygon 6               42  5.58317e+04      7.46e-05
polygon 7               67  1.31354e+06      1.75e-03
polygon 8               15  4.46420e+03      5.96e-06
polygon 9               14  5.46674e+03      7.30e-06
polygon 10              37  5.26194e+03      7.03e-06
polygon 11              53  3.44003e+04      4.59e-05
polygon 12              74  5.82234e+04      7.78e-05
polygon 13              69  5.63134e+04      7.52e-05
polygon 14             143  1.45139e+05      1.94e-04
polygon 15             165  3.38736e+05      4.52e-04
polygon 16             130  9.40465e+04      1.26e-04
polygon 17              19  1.80977e+03      2.42e-06
polygon 18              16  2.01046e+03      2.69e-06
polygon 19              93  4.30642e+05      5.75e-04
polygon 20              90  4.15092e+05      5.54e-04
polygon 21             721  1.92795e+06      2.57e-03
polygon 22             330  1.11896e+06      1.49e-03
polygon 23             115  9.28394e+05      1.24e-03
polygon 24              37  1.01705e+04      1.36e-05
polygon 25              25  1.66227e+04      2.22e-05
polygon 26              10  2.14507e+03      2.86e-06
polygon 27             190  2.02489e+05      2.70e-04
polygon 28             175  9.25904e+05      1.24e-03
polygon 29            1990  9.99217e+06      1.33e-02
polygon 30 (hole)        3 -1.06765e+00     -1.43e-09
polygon 31              38  2.42492e+04      3.24e-05
polygon 32              24  6.35239e+03      8.48e-06
polygon 33              53  6.35791e+05      8.49e-04
polygon 34              41  1.60161e+04      2.14e-05
polygon 35              22  2.54368e+03      3.40e-06
polygon 36              30  1.08382e+04      1.45e-05
polygon 37             327  2.16921e+06      2.90e-03
polygon 38             111  6.62927e+05      8.85e-04
polygon 39              90  1.15991e+05      1.55e-04
polygon 40              98  6.26829e+04      8.37e-05
polygon 41             415  3.25384e+06      4.35e-03
polygon 42             222  1.51142e+06      2.02e-03
polygon 43             107  6.33039e+05      8.45e-04
polygon 44               7  2.48299e+03      3.32e-06
polygon 45              17  3.28303e+04      4.38e-05
polygon 46              26  8.34758e+03      1.11e-05
polygon 47             177  4.67446e+05      6.24e-04
polygon 48              16  3.19460e+03      4.27e-06
polygon 49              15  4.87296e+03      6.51e-06
polygon 50              66  1.61841e+04      2.16e-05
polygon 51             149  5.63430e+06      7.53e-03
polygon 52             609  2.62570e+07      3.51e-02
polygon 53               8  7.82256e+03      1.04e-05
polygon 54             976  2.33447e+07      3.12e-02
polygon 55              55  8.25379e+04      1.10e-04
polygon 56             976  2.33447e+07      3.12e-02
polygon 57              61  3.33449e+05      4.45e-04
polygon 58               6  1.68410e+04      2.25e-05
polygon 59               4  9.45963e+03      1.26e-05
polygon 60              46  6.99702e+05      9.35e-04
polygon 61              13  7.00873e+04      9.36e-05
enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
                     (53380 x 33890 units)
Window area = 748741000 square units
Fraction of frame area: 0.414

1st Order SPPA

rescale to convert from m to km

airbnbSG_ppp.km <- rescale(airbnbSG_ppp, 1000, "km")

Kernel Density Estimation (automatic bandwidth)

kde_airbnbSG.bw <- density(airbnbSG_ppp.km,
                              sigma=bw.diggle,
                              edge=TRUE,
                            kernel="gaussian") 

other methods to do: bw.CvL(), bw.scott(), bw.ppl() other smoothing kernels: epanechnikov, quartic, disc

plot(kde_airbnbSG.bw)

KDE (fixed bandwidth)

kde_airbnbSG_600 <- density(airbnbSG_ppp.km, sigma=0.6, edge=TRUE, kernel="gaussian") #0.6 sigma = 600 m for bandwidth. 600m = 0.6km
plot(kde_airbnbSG_600)

sigma: range from 0.1 to 10 different kernels

KDE (adaptive bandwidth)

kde_airbnbSG_adaptive <- adaptive.density(airbnbSG_ppp.km, method="kernel")
plot(kde_airbnbSG_adaptive)

convert KDE output into grid (do for all)

gridded_kde_airbnbSG_bw <- as.SpatialGridDataFrame.im(kde_airbnbSG.bw)

# convert grid output into raster
kde_airbnbSG_bw_raster <- raster(gridded_kde_airbnbSG_bw)

# assign projection systems
projection(kde_airbnbSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_airbnbSG_bw_raster

class      : RasterLayer 
dimensions : 128, 128, 16384  (nrow, ncol, ncell)
resolution : 0.4170614, 0.2647348  (x, y)
extent     : 2.663926, 56.04779, 16.35798, 50.24403  (xmin, xmax, ymin, ymax)
crs        : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs 
source     : memory
names      : v 
values     : -1.331664e-13, 633.9967  (min, max)

final visualisation tmap

tm_shape(kde_airbnbSG_bw_raster) + 
  tm_raster("v") +
  tm_layout(legend.position = c("right", "bottom"), frame = FALSE)

2nd Order SPPA G/F/K/L Function

G Function estimation

G_airbnb = Gest(airbnbSG_ppp, correction =c("rs", "Hanisch"))
plot(G_airbnb, xlim=c(0,500))

other correction: “best” Gest() - corrections: “none”, “rs”, “km”, “Hanisch”, “best”, “all” Fest() - corrections: “none”, “rs”, “km”, “cs”, “best”, “all” Kest() - corrections: “none”, “border”, “bord.modif”, “isotropic”, “Ripley”, “translate”, “translation”, “rigid”, “none”, “good”, “best”, “all” Lest() - same as Kest()

Complete Spatial Randomness Test

G_airbnb.csr <- envelope(airbnbSG_ppp, Gest, nsim = 999)
# G_airbnb.csr <- envelope(airbnbSG_ppp, Gest, correction="all", nsim = 999)

nsim: number btwn 0 to 1000

plot to show

plot(G_airbnb.csr)

visualisation

tmap_mode("plot")
tm_shape(mpsz_sf) + 
  tm_borders(alpha = 0.5) +
tm_shape(airbnb_sf) + 
  tm_dots(col = "red", size = 0.05) +
  tm_layout(main.title = "Airbnbs",
            main.title.position = 'center',
            main.title.size = 0.8,
            frame = TRUE)

occurrence <- airbnb_sf %>% count(airbnb_sf$neighbourhood)

occurrence <- occurrence[order(-occurrence$n),]

occurrence

Simple feature collection with 43 features and 2 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 7406.989 ymin: 25352.92 xmax: 44064 ymax: 48321.84
Projected CRS: SVY21 / Singapore TM
# A tibble: 43 × 3
   `airbnb_sf$neighbourhood`     n                                      geometry
   <chr>                     <int>                              <MULTIPOINT [m]>
 1 Kallang                     353 ((29508.14 32976.25), (29791.93 33174.18), (…
 2 Downtown Core               325 ((28945.03 28564.31), (28952.82 28557.68), (…
 3 Outram                      255 ((28290.63 29757.42), (28596.68 29346.08), (…
 4 Rochor                      214 ((29383.5 31768.77), (29521.51 31404.98), (2…
 5 Novena                      182 ((28059.14 33734.79), (28062.48 33523.59), (…
 6 Queenstown                  176 ((20006.09 30435.35), (20082.89 30453.04), (…
 7 Bukit Merah                 163 ((24405.42 28144.15), (24581.26 28085.54), (…
 8 River Valley                157 ((27291.23 31642.72), (27473.75 31495.65), (…
 9 Geylang                     145 ((32942.57 32745.19), (33084.45 32747.7), (3…
10 Bedok                       120 ((36040.87 33769.19), (36077.61 33347.9), (3…
# … with 33 more rows

kallang has the most airbnbs

K cross function

NetSPPA

kallang airbnbs:

kallang_airbnb_sf <- filter(airbnb_sf, neighbourhood == "Kallang")
kallang_airbnb_sf

Simple feature collection with 353 features and 4 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 29508.14 ymin: 31119.73 xmax: 33791.74 ymax: 34536.48
Projected CRS: SVY21 / Singapore TM
# A tibble: 353 × 5
   name                          neigh…¹ room_…² price            geometry
 * <chr>                         <chr>   <chr>   <dbl>         <POINT [m]>
 1 EnsuiteA w Balcony view/acce… Kallang Privat…    57 (31346.67 32390.22)
 2 Ensuite M1B within Heritage … Kallang Privat…    50    (30798 32708.67)
 3 EnsuiteK w Balcony view/acce… Kallang Privat…    57 (31239.83 32106.04)
 4 Heritage apartment: Master r… Kallang Privat…    69 (30505.31 32834.72)
 5 Lavender QueenBed *rmFr: 5m … Kallang Privat…    59 (31368.93 32169.07)
 6 Life Impact Coaching          Kallang Entire…   130 (32498.52 32598.12)
 7 1-Pax Small Private Room, Sh… Kallang Privat…    78 (31000.55 32666.66)
 8 2-pax Small Private Room, Sh… Kallang Privat…   110 (31000.55 32666.66)
 9 1-Pax Small Private Room, Sh… Kallang Privat…    78 (31000.55 32666.66)
10 1950s cozy window singlebed … Kallang Privat…    49 (31207.55 32608.05)
# … with 343 more rows, and abbreviated variable names ¹​neighbourhood,
#   ²​room_type

plot(network_sf)

kallang airbnb plot

tmap_mode('plot')
tm_shape(kallang_airbnb_sf) + 
  tm_dots() + 
  tm_shape(network_sf) +
  tm_lines()

all airbnbs plot

tmap_mode('view')
tm_shape(airbnb_sf) + 
  tm_dots() + 
  tm_shape(network_sf) +
  tm_lines()

tmap_mode('plot')

convert to single linestring

network_linestring <- network_sf %>% st_cast("LINESTRING")

save to rds

network_linestring_rds <- write_rds(network_linestring, "data/network_linestring.rds")

preparing lixels objects

lixels <- lixelize_lines(network_linestring, 
                         700, 
                         mindist = 350)

save to rds

lixels_rds <- write_rds(lixels, "data/lixels_rds.rds")

read the rds file

lixels_rds <- read_rds("data/lixels_rds.rds")

generate line centre points

samples <- lines_center(lixels_rds)

write to rds

samples_rds <- write_rds(samples, "data/samples_rds.rds")

read the rds file

samples_rds <- read_rds("data/samples_rds.rds")

perform NetKDE

densities <- nkde(network_linestring, 
                  events = kallang_airbnb_sf,
                  w = rep(1,nrow(kallang_airbnb_sf)),
                  samples = samples_rds,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)

write to rds

densities_rds <- write_rds(densities, "data/densities_rds.rds")

read the rds file

densities_rds <- read_rds("data/densities_rds.rds")

samples_rds$density <- densities_rds
lixels_rds$density <- densities_rds

# rescaling to help the mapping
samples_rds$density <- samples_rds$density*1000
lixels_rds$density <- lixels_rds$density*1000

tmap_mode('view')
tm_shape(lixels_rds)+
  tm_lines(col="density")+
tm_shape(kallang_airbnb_sf)+
  tm_dots()

kfun_kallang_airbnb <- kfunctions(network_linestring, 
                             kallang_airbnb_sf,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 50, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05,
                             agg = 40)

kfun_kallang_airbnb$plotk

for all airbnbs

perform NetKDE for all airbnbs

all_densities <- nkde(network_linestring, 
                  events = airbnb_sf,
                  w = rep(1,nrow(airbnb_sf)),
                  samples = samples_rds,
                  kernel_name = "quartic",
                  bw = 300, 
                  div= "bw", 
                  method = "simple", 
                  digits = 1, 
                  tol = 1,
                  grid_shape = c(1,1), 
                  max_depth = 8,
                  agg = 5, #we aggregate events within a 5m radius (faster calculation)
                  sparse = TRUE,
                  verbose = FALSE)

write to rds

all_densities_rds <- write_rds(all_densities, "data/all_densities_rds.rds")

read rds

all_densities_rds <- read_rds("data/all_densities_rds.rds")

samples_rds$density <- all_densities_rds
lixels_rds$density <- all_densities_rds

# rescaling to help the mapping
samples_rds$density <- samples_rds$density*1000
lixels_rds$density <- lixels_rds$density*1000

tmap_mode('view')
tm_shape(lixels_rds)+
  tm_lines(col="density")+
tm_shape(airbnb_sf)+
  tm_dots()

kfun_all_airbnbs <- kfunctions(network_linestring, 
                             airbnb_sf,
                             start = 0, 
                             end = 1000, 
                             step = 50, 
                             width = 50, 
                             nsim = 50, 
                             resolution = 50,
                             verbose = FALSE, 
                             conf_int = 0.05,
                             agg = 1500)

kfun_all_airbnbs$plotk

network cross K function spNetwork

crossk_kallang <- cross_kfunctions(network_linestring, kallang_airbnb_sf, mrt_sf, start = 0, end = 5000, step = 50, width = 1000, nsim = 50, verbose = FALSE, agg = 100)

crossk_kallang$plotk

crossk_kallang <- cross_kfunctions(network_linestring, mrt_sf, kallang_airbnb_sf, start = 0, end = 5000, step = 50, width = 1000, nsim = 50, verbose = FALSE, agg = 100)

crossk_kallang$plotk

kcross

Kcross(airbnbSG_ppp.km, )